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1.0 Introduction 

Let Xj,.. ,X N be a sample of N observations taken from a bivariate normal distribution with unknown mean 
vector and covariance matrix E, and let X and S be the respective 2x1 sample mean vector and 2x2 sample 
covariance matrix calculated from the X-. Given a desired containment probability 0 and a level of confidence y , 
the problem addressed is to find a region about X that contains at least 1000% of the X-distribution with probability 

7- 

When S~ ] exists, an ellipse R about X for any positive number c may be defined by 

R = R(X,S t c) = {x: * c} (1-1) 

For a future observation X distributed according to N 2 (pX)> the probability that X t R is given by 


He) = 


r/« 


{x-X )'^ 1 (x-Jf) 


dx 


( 1 - 2 ) 


For each new sample of N observations, R , which depends on X and S, is a random region in 2-dimensional 
Euclidean space; thus for a fixed general c, 1(c) is a random variable taking values in (0,1). In particular, we seek 
c* (not depending on the unknown (x or L), such that 


P{ J(c+) ;> p} = y (1-3) 

The corresponding ellipse R(X ,S,c*), known as a tolerance region , solves the problem. 

Based on the work of John (1963) for a general /7-variate normal distribution, the following approximation 
for c* is given by Chew (1966): 


(M) 

VL^iiN-Dp) 

where u $(p,\) is the ^-percentage point of the noncentral chi-squared distribution 1 with p degrees of freedom and 
noncentrality parameter X, and u 2 (m) is the (1 -y ) -percentage point of the central chi-squared distribution with m 
degrees of freedom. Equation (1-4) is easier to evaluate than more precise but complicated expressions, such as 
given by Siotani (1964). Chew states, "the approximation is good if 1/N 2 is negligible"; however, it appears (see 
section 3) that for the bivariate normal case (p = 2), c underestimates c * by a factor 1 - A/N where A depends on 
0 and 7. 

For general values of p y (1-4) has stood the test of time; for example it was cited and used by Rode and 
Chincilli (1988) in their paper on transforming clinical laboratory measurements. When p — 2, however, it is 
feasible to significantly improve the approximation by direct calculation of 1(c) within a Monte-Carlo simulation of 
values of X and S (see section 2). Estimation of A by comparing the resulting more accurate estimates of c with 
c makes it feasible to use a corrected form of (1-4) to obtain accurate easily computed tolerance regions. 


1 Chew defines the noncentrality parameter "in accordance with that in Wilks" (1962); i.e. a noncentral chi- 

squared random variable with m degrees of freedom and noncentrality parameter X is distributed as Z 2 -f F, where 
Z — N(X 7/2 ,I) and Fhas a central chi-squared distribution with m-1 degrees of freedom. 
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2,0 Monte-Carlo Estimation of c* 


For X - N(fi t E) and any level of confidence, it can be shown that as //becomes large, c approaches c 0 
= -2log(l-0). This is because c 0 satisfies 


P^X-ijO'E-MX-h) *c 0 } = 1-e 2 =P 


( 2 - 1 ) 


(e.g., see Cramer (1963)) and X and S converge in probability to ^ and E as N increases. For finite N, the solution 
to (1-3) is c* = Kc 0 for some K > 1. 

Let L 1/2 be a "square-root" of E in the sense that Z 1/2 (L 1/2 ) t = E. By making the transformation y = 
E ~ 1/2 ( x-fx) in (1-2), it can be shown that the solution to (1-3) is the same as when fi = 0, E = / and X and S are 
obtained from a sample of N observations from the N(0,7)-distribution. As a result, it will be henceforth assumed 
that fx - 0 and E = I. 

For each combination of N = 10, 40(5), 50, 0 =.90, .95, .99, .999 and selected values of c in the range 
c - Kc 0 (1 < K £ 7.5), 1000 realizations of X and S were randomly generated taking ^ = 0 and E= I. (This 
can be done without generation of the individual observations; see Odell and Feiveson (1966)). For each X and S , 
1(c) was then calculated by numerical integration (see appendix). 

With 0 fixed, Q(c) = P{I(c) ^ 0} is a monotonic increasing function of c, with c* being the root Q(c ) 
= y. From the simulation, for each trial value of c , say c it the observed proportion of times, q i9 that l(c0 exceeds 
0, is an estimate of Q(c0. For each N and 0 , an interpolating quadratic function was fitted to the points (y if cj, 
where y t = -logfl-qj; (.80 < q i <> .999), then set equal to y y = -log(l-y) to solve for c y , the estimate of for 
y = .90, .95 and .99. As an example, a plot of y { vs c i along with the interpolating quadratic function is shown for 
N = 10 and 0 = .99 in figure 1. The three horizontal lines represent the values of y y which define c y . 



Figure 1. y t = -log(l-q0 c- for N = 10; 0 = .99 
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Inin Him it 



Originally, X and S were kept fixed as c was varied for given N and 0\ however, it was noticed that unnatural 
patterns in plots of y vs c would often result. Consequently, it was decided to avoid all dependence between results 
by regenerating X and S for each value of c , despite the extra computational effort. 

3.0 Accuracy of Point Estimates 

The major sources of error in c y are (1) the numerical integration used to compute 1(c) and (2) the process 
of fitting an interpolating polynomial to the q- and solving for c. A check against the first error was made by using 
a sufficiently small step size so that values of 1(c) resulted in 1 (to within 5 decimal places) as c was made 
arbitrarily large. The second error contains a random component induced by the binomial distribution of the q i and 
some bias due to inverting the estimated y vs c relationship as well as possible model error (the true relationship 
may not be quadratic). Because the obtained fits were so tight (e.g., see fig. 1), the bias in c y was considered 
negligible. 

To estimate the variance of the random error, a replicate of the entire simulation was made. Under the 
assumption of constant coefficient of variation, differences between the results were used to estimate a CV of about 
1.5% for individual values of c y . For the final smoothing described below, results of the two runs were averaged, 
further reducing the error CV by a factor of ^ 2. 

4.0 Final Results 

Values of c c 95 , and c 99 for each Vand 0 are shown in table 1 along with corresponding values of c 
obtained from (1-4). A comparison reveals that the latter tend to be smaller by a factor of 1 - A/N where A depends 
on 0 and 7 , thus suggesting smoothing the c y using c as a concomitant variable. Given c , one may then better 
estimate c* by 

c f = c[N/(N-A) ]. (4-1) 

Using the data in table 1 , estimated values of A for various 0 and 7 were obtained by regression through the origin 
of 1 - c y /c y against 1/N. These values are shown in table 2. Although there were only eight values of N for each 
of the nine regressions, the fits were almost exact. Standard errors of ^-estimates ranged from .04 to .22, 
corresponding to pertubations in c' between 0.4 and 2.3 percent for N = 10 , and between 0.08 and 0.45 percent 
for N = 50. By contrast, errors in uncorrected c are about 50% for N = 10 and 5 - 10% for N — 50. 


Table 1. Values of c y and c 
7 -Confidence Tolerance Ellipsoid of Content 0\ 
Bivariate Normal Distribution 


N 

P 

11 

.90 

d y 

II 

.95 

S y 

ll 

.99 

10 

.900 

12.53 

8.38 

15.45 

9.69 

24.93 

12.97 

10 

.950 

17.07 

10.89 

21.40 

12.60 

34.69 

16.86 

10 

.990 

28.05 

16.72 

35.46 

19.34 

58.85 

25.89 

10 

.999 

45.42 

25.06 

57.25 

28.99 

91.78 

38.81 
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Table 1 (Cont.), Values of c y and c 
y-Confidence Tolerance Ellipsoid of Content 0 ; 
Bivariate Normal Distribution 




Y = 

.90 

Y = 

.95 

Y = 

.99 

N 

p 

S 

S y 

S 

s 

S 

S y 

15 

.900 

9.23 

7.26 

10.68 

8.12 

14.33 

10.13 

15 

.950 

12.47 

9.43 

14.67 

10.56 

21.01 

13.17 

15 

.990 

20.23 

14.50 

23.65 

16.22 

32.80 

20.24 

15 

.999 

31.75 

21.71 

37.90 

24.29 

52.91 

30.31 

20 

.900 

7.91 

6.72 

8.93 

7.38 

11.36 

8.87 

20 

.950 

10.44 

8.74 

11.78 

9.60 

15.09 

11.54 

20 

.990 

16.86 

13.43 

19.20 

14.75 

25.43 

17.74 

20 

.999 

26.57 

20.15 

30.21 

22.14 

40.58 

26.63 

25 

.900 

7.25 

6.39 

8.00 

6.94 

9.95 

8.16 

25 

.950 

9.66 

8.32 

10.68 

9.03 

13.07 

10.61 

25 

.990 

15.31 

12.78 

17.09 

13.88 

21.57 

16.30 

25 

.999 

24.19 

19.19 

27.09 

20.85 

34.43 

24.49 

30 

.900 

6.83 

6.17 

7.45 

6.65 

8.97 

7.68 

30 

.950 

9.04 

8.03 

9.89 

8.65 

11.91 

9.99 

30 

.990 

14.23 

12.34 

15.78 

13.30 

19.65 

15.36 

30 

.999 

22.26 

18.49 

24.58 

19.92 

30.19 

23.01 

35 

.900 

6.53 

6.01 

7.01 

6.44 

8.27 

7.35 

35 

.950 

8.58 

7.82 

9.32 

8.38 

11.11 

9.56 

35 

.990 

13.52 

12.02 

14.65 

12.87 

17.50 

14.69 

35 

.999 

21.33 

18.02 

23.33 

19.29 

27.52 

22.01 

40 

.900 

6.30 

5.89 

6.78 

6.28 

7.91 

7.09 

40 

,950 

8.36 

7.66 

8.98 

8.16 

10.40 

9.23 

40 

.990 

13.14 

11.78 

14.20 

12.55 

16.78 

14.18 

40 

.999 

20.31 

17.67 

21.97 

18.83 

26.04 

21.27 

50 

.900 

6.02 

5.71 

6.42 

6.04 

7.33 

6.73 

50 

.950 

7.90 

7.43 

8.40 

7.86 

9.70 

8.75 

50 

.990 

12.52 

11.43 

13.36 

12.09 

15.40 

13.46 

50 

.999 

19.03 

17.15 

20.41 

18.13 

23.75 

20.19 


Table 2. Values of A for Correcting 
Non-central Chi-Squared Approximation for 
Bivariate Normal Tolerance Regions. 

Y 

0.90 0.95 0.99 

0.900 3.153 3.543 4.553 

P 0.950 3.521 3.994 5.103 

0.990 4.093 4.606 5.800 

0.999 4.725 5.254 6.334 

Correction is c f = c [N/(N -A)] . 

As an example, for N = 10, a 90%-tolerance region (7 = .90) that contains at least 99% of the population (0 
= .99) is found by first computing the chi-squared approximation (1-4), giving c = 16. 72 , and then correcting 
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it with equation (4-1). Table 2 gives A = 4.093; hence c* = 16. 72[10/(10 - 4.093)] = 28.31. The desired tolerance 
region is the ellipse {x:( x-X)'S~ J (x-X) ^ 28.31}, 

5.0 Concluding Remarks 

This paper has illustated how Monte-Carlo simulation, along with simple regression modelling, can be used 
to improve a theoretical approximation for a useful special case. The approximation is easily obtained if one has 
access (through software or tables) to the percentage points of the central and non-central chi-squared distributions. 
Correction to more accurate values for bivariate normal tolerance regions is readily accomplished for conventional 
values of /? and 7 using the appropriate value of A in table 2. 

If one does not have a ready means of obtaining non-central chi-squared percentage points u'p(p,\), an 
approximation given in Abramowitz and Stegun (1966) provides even greater simplification of computation with little 
loss of accuracy when p = 2. Abramowitz and Stegun give u^(p,\) « (1 +b)ug(p*) where p* = a/(l + b), a = 
p + X and b — \/(p + \). Here, p = 2 and X = 2/N, hence 1 + b — (N+2)/(N+ 1) and p* — 
2(N+l) 2 /[N(N+2)]~ 2 so that for larger values of N, one may simply use 

up (p, X) * 


(5-1) 
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Appendix 


Obtaining 1(c) by Numerical Integration 
Equation (1-2) may be rewritten 

*2H sr*u 2 > 

I(c)=j f(x 2 ) j f(x 1 \x 2 )dx 1 dx 2 (A-l) 

X 2L 


where f( x 2) is the density of x 2 ; i.e. N(0,1) and/fx ; \x^) is the density of the conditional distribution of x l given x 2 , 
which is also standard normal, since E= 7. The limits x 2H and x 2L are given by X 2 ±(cS 22 ) 1/2 and for fixed x 2 , the 
limits of Xj are given by 




S 2l ( ~ X 2 ) ±yJ\S\ [cS 2 2 (-^2 ^2 ) ^ 


(A-2) 


where S = (S- ). The inner integral in (A-l) is easily computed as w bere $ is the standard 

normal cumulative distribution function for which good approximations are available. 
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